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ABSTRACT 

A  seven  level  primitive  equation  equatorial  ocean  model 
with  temperature  and  salinity  is  integrated  over  a  30  year 
period  to  a  quasi-steady  state.   Placement  of  four  of  the 
seven  levels  of  the  model  in  the  first  100  meters  gives 
sufficient  vertical  resolution  to  resolve  most  major  features 
of  equatorial  oceans.   The  model  equivalents  of  the  South 
Equatorial  Current,  the  Cromwell  Current,  vertical  upwelling 
in  the  equatorial  area,  Ekman  divergence  of  mass  toward  each 
pole,  and  a  thermocline  were  all  simulated. 

A  major  conclusion  was  salinity  probably  has  little 
effect  on  the  general  circulation  patterns  in  the  equatorial 
region.   Additionally,  the  use  of  a  vertical  eddy  diffusion 
coefficient  for  momentum  that  is  dependent  upon  vertical 
stability,  vertical  current  shear,  and  depth  will  probably 
be  required  to  accurately  model  the  vertical  structure  of 
the  currents. 
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I.   INTRODUCTION 

Numerical  baroclinic  ocean  models  are  advancing  to  the 
state  where  oceanic  forecasting  may  become  feasible.   The 
motivation  for  the  development  of  sophisticated  ocean  models 
stems  largely  from  the  need  to  accurately  specify  the  lower 
boundary  conditions  for  baroclinic  atmospheric  models  used 
in  long-range  weather  prediction.   The  importance  of  the 
ocean  in  determining  global  climate  is  underscored  in  a 
recent  pioneering  work  by  Manabe  (1969)  and  Bryan  (1969)  in 
which  numerical  ocean  and  atmospheric  models  were  coupled 
together  allowing  interaction  between  the  two.   By  comparing 
the  results  of  the  ocean-air  model  with  those  from  an  atmos- 
pheric model  with  a  static  ocean,  they  found  the  ocean  plays 
an  important  role  in  the  poleward  transport  of  heat  and 
helps  determine  the  thermal  structure  and  vertical  wind  shear 
of  the  troposphere.   Additionally,  they  found  that  the  ocean 
affects  the  hydrology  of  the  tropics  and  influences  the 
strength  of  the  Hadley  circulation.   The  important  role  of 
the  ocean  in  influencing  the  transient  behavior  of  the 
atmosphere  has  been  recently  shown  by  Lambertson  (1972)  in 
a  coupled  ocean-atmospheric  model.   Because  of  the  demon- 
strated importance  of  the  ocean  in  determining  the  climate 
and  the  variability  of  the  atmosphere,  and  since  atmospheric 
models  are  more  highly  developed  than  their  ocean  counter- 
parts, it  is  clear  that  future  improvement  in  long-range 
weather  prediction  is  closely  linked  with  the  development  of 
better  ocean  models. 


Early  oceanic  prediction  models  were  one-dimensional  in 
the  vertical  with  the  net  heat  flux  and,  in  some  cases,  the 
wind  velocity  or  stress  being  specified  at  the  surface. 
The  thermal  structure  was  assumed  to  be  known  initially 
and  responded  to  changes  in  the  rates  of  surface  heating. 
The  depth  of  the  surface  layer  was  determined  by  a  mechan- 
ical energy  balance,  and  integration  of  the  thermodynamic 
equation  gave  the  thermal  profile  at  a  later  time.   Although 
adequate  for  very  short  range  purposes,  these  one-dimensional 
models  are  not  satisfactory  for  long  range  predictions 
because  no  mechanisms  for  horizontal  advection  or  diffusion 
of  heat  or  momentum  are  present.   Any  location  where  the 
climatological  yearly  heat  flux  is  not  zero  will  display 
either  a  continual  warming  or  cooling  trend.   Tropical 
areas  would  continue  to  warm  and  middle  latitude  areas 
would  continue  to  cool.   For  long  term  predictions,  there- 
fore, ocean  models  based  on  heat  balances  alone  will  not 
be  satisfactory. 

Although  Bryan's  dynamical  ocean  model  is  adequate  for 
the  purpose  of  global  climate  investigation,  it  is  believed 
that  the  horizontal  grid  spacing  of  5   longitude  and  5 
latitude  and  a  vertical  spacing  where  the  first  level  is 
at  67  meters  is  too  coarse  for  adequate  oceanic  forecasting 
purposes.   One  strong  point  of  the  present  model  is  felt 
to  be  its  superior  horizontal  resolution  of  2.5   latitude 
and  3-0   longitude  and  good  vertical  resolution  where  four 
of  its  seven  levels  are  in  the  first  100  meters. 
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Because  of  the  importance  of  equatorial  air-sea  inter- 
action (Bjerknes,  1966)  and  the  interesting  oceanographic 
phenomena  that  occur  in  the  equatorial  region,  it  was 
decided  to  model  the  equatorial  domain.   The  large  vertical 
current  shears  and  sharp  thermocline  which  exist  in  the 
tropics  should  be  simulated  by  an  ocean  model  that  is  to 
be  used  with  an  atmospheric  forecasting  model. 
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II.   DESIGN  OF  THE  MODEL 

A.   MATHEMATICAL  STATEMENT  OP  PROBLEM 

The  idealized  domain  is  60  wide,  extends  from  21.25  S 
to  21.25  N,  and  has  a  uniform  depth  of  3000  meters.   Place- 
ment of  the  lateral  walls  at  approximately  20   latitude 
in  both  hemispheres  satisfies  two  important  requirements: 
first,  the  curl  of  the  observed  wind  stress  is  almost  zero 
at  these  locations;  and  second,  the  area  averaged  flux  of 
salt  as  determined  from  climatological  evaporation  and 
precipitation  data  is  approximately  zero  when  averaged  from 
21.25  S  to  21.25  N.   Seven  vertical  levels  placed  at  5,  20, 
50,  100,  250,  700,  and  2000  meters  are  used.   The  upper 
levels  are  more  closely  spaced  to  provide  better  resolution 
in  the  near  surface  domain  where  the  largest  variation  in 
dependent  variables  is  found.   Horizontally  the  grid  points 
are  spaced  3.0  degrees  longitude  and  2.5  degrees  latitude. 
The  horizontal  and  vertical  domains  are  pictured  in  figures 
2.1  and  2.2  respectively. 

This  model  is  similar  to  the  primitive  equation  models 
of  Haney  (1971),  Bryan  (1969),  and  Bryan  and  Cox  (1968) 
except  for  the  introduction  of  salinity  as  an  explicit  de- 
pendent variable.   The  solution  of  the  equations  of  motion 
involves  the  use  of  a  streamfunction  for  the  vertical  mean 
motion.   This  somewhat  artificial  constraint  allowed  the  use 
of  a  large  time  step  and  saved  considerable  computation  time 


21.25  N 


21.25  S 


Figure  2.1   Horizontal  domain  of  the  model.   Horizontal 
velocity  components  u  and  \r  were  carried  between  grid 
points.   The  variables  T,  S,  f,  ¥,w  were  carried  on  the 
grid  points.   Grid  spacing  was  3-0  degrees  longitude  and 
2.5  degrees  latitude. 
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Figure  2.2   Vertical  structure  of  the  model  showing  seven 
levels  with  dense  packing  near  the  surface.   All  variables 
were  stored  at  these  levels  except  vr;  the  vertical  velocity, 
which  was  found  at  intermediate  levels  (not  shown,  see 
figure  2.8). 
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The  model  is  based  upon  the  finite  difference  forms  of 
the  primitive  equations  which  are  integrated  in  time  from 
a  known  set  of  initial  conditions.   The  governing  equations 
are  the  two  horizontal  equations  of  motion,  the  hydrostatic 
equation,  the  continuity  equation,  thermodynamic  equa- 
tion, the  salinity  equation,  and  the  equation  of  state 


du   _  .       » 


l£  -v  (2.H+  m    Wshhj  ♦  FX)  (2.1) 


di  p0a  c©&«f 


<tU  Ah71T  +  feH^T  ^    CtCT^V 
dt  t)2*  ' 


1 


where  the  symbols  are  as  follows: 
*fc   time 
X    longitude 

latitude 
t£    height 
Jft,   angular  speed  of  the  earth 
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(2.2) 


(2.3) 


4*  } 


(2.5) 


iSm    A*V2S    +    k$££    +    CslT^S).  (2.6) 

p»   ?o[»-c*(7-T6)  +   ^S-So)"],  (2.7) 


0.        radius  of  the  earth 

q        acceleration  of  gravity 

u        zonal  current  component 

\T        meridional  current  component 

ty        horizontal  velocity 

W*        vertical  speed 

T        temperature 

"To       constant  reference  temperature 

'S        salinity 

S&       constant  reference  salinity 

Ct(Tj^)    rate  of  change  of  T  due  to  instantaneous 

convective  adjustment 
C$OjS)    rate  of  change  of  &  due  to  instantaneous 

convective  adjustment 
p        density 
^o        density  of  water  at  reference  salinity  $0  and 

temperature  To 
«l        coefficient  of  thermal  expansion  of  sea  water 
&        coefficient  of  saline  contraction  of  sea  water 
fcn^e     vertical  eddy  diffusion  coefficients'  for  heat 

and  salt,  respectively 
Ah,&$    lateral  eddy  diffusion  coefficients  for  heat  and 

salt,  respectively 
rh  eastward  frictional  force  per  unit  volume 

P*J       northward  frictional  force  per  unit  volume 
V        horizontal  gradient  operator. 


16 


The  model  circulation  is  determined  primarily  by  the 
surface  wind  stress,  and  to  a  lesser  extent  by  the  fluxes  of 
heat  and  salt  at  the  surface.  Figure  2.3  shows  the  annual 
latitudinal  mean  zonal  stress  for  all  oceans  combined 
(Hellerman,  1967)  which  was  used  in  this  model.   The  distri- 
bution is  practically  symmetric  except  for  the  occurrence  of 
a  single  minimum  of  the  easterly  stress  north  of  the  equator, 

The  net  heat  flux  at  the  upper  boundary  is  calculated 
according  to  Haney  (1971)  from  the  equation 


Qs  QL4  Qi.tTk-Ti), 


(2.8) 


where : 

Q       is  the  net  heat  flux  into  the  ocean, 

Qj,  is  that  component  of  heat  flux  that  would  exist 
if  there  were  no  air-sea  temperature  difference 
(this  is  primarily  radiative  fluxes), 

Q^.      is  the  contribution  due  to  the  existence  of  a 
temperature  differential  between  the  ocean  and 
the  atmosphere  (this  is  primarily  sensible  and 
latent  heat  fluxes),  and 

TK      is  a  prescribed  surface  air  temperature. 

Haney' s  results,  which  were  based  on  climatological  data 
after  London  (1957),  were  averaged  across  the  equator  to 
obtain  the  symmetric  patterns  of  Qj_,  Q^j  and  ~^k  shown  in 
figure  2.k.  This  form  of  boundary  condition  essentially 
couples  the  ocean  surface  temperature  to  a  prescribed 
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Figure  2.3   The  annual  latitudinal  mean  zonal  stress  for 
all  oceans  combined  after  Hellerman  (1967). 


18 


UJ 
o 

ZD 
H 

h 
< 


LYDAY"1    or     LYDAY^'dEG"1 
0  100  200 


Figure  2.4  The  atmospheric  surface  air  temperature 
(bottom  scale  in  degrees  Celsius)  and  the  heat  flux 
components  Q\_   and  Qj_  in  lyday"!  and  lyday~ldeg"*l  (top  scale) 
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atmospheric  temperature.   To  accomplish  a  net  heat  flux  of 
zero  requires  that  T^  be  less  than  "T^  over  at  least  some 
portions  of  the  domain. 

The  boundary  condition  on  salinity  is  quite  different 
from  that  on  temperature  because  the  downward  flux  of 
salinity  is  determined  by  the  rate  of  removal  of  fresh  water 
from  the  surface  of  the  ocean.   Rather  than  specifying  the 
salinity  at  the  surface,  the  downward  flux  of  salt  was  speci- 
fied from  the  zonal  mean  climatological  rates  of  precipita- 
tion and  evaporation  after  Budyko  (1963).   These  curves  are 
shown  in  figure  2.5.   The  data  were  averaged  across  the 
equator  to  obtain  the  symmetric  curve  shown  by  the  dotted 
line.   In  this  way,  the  only  asymmetric  boundary  condition 
was  the  surface  wind  stress.   This  curve  has  the  property 
that  the  meridionally  integrated  net  rate  of  evaporation 
minus  precipitation  is  zero  and,  since  no  flux  of  salt  was 
allowed  at  the  lateral  boundaries  or  ocean  floor,  was  suffi- 
cient to  guarantee  that  the  volume  averaged  salinity  of  the 
ocean  remained  constant  in  time.   This  does  not  imply,  how- 
ever, that  no  variations  in  surface  salinity  patterns  would 
tend  to  occur.   Clearly,  near  the  equator  where  E-P  has  a 
local  maximum,  a  local  maximum  of  salt  would  tend  to  be 
formed.   In  the  region  of  the  model  where  precipitation 
exceeded  evaporation  the  most,  a  local  minimum  would  tend  to 
be  formed.   Whether  the  horizontal  and  vertical  gradients  of 
salt  caused  by  the  upper  boundary  condition  would  modify  the 
model  circulation  was  an  important  question  to  be  answered  by 
this  investigation. 
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Figure  2.5   Evaporation  and  precipitation  rates  (solid 
lines)  after  Budyko  (1963).   The  net  rate  of  evaporation 
minus  precipitation  averaged  across  the  equator  is  shown 
by  the  dashed  line. 
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Details  of  the  space  differencing  of  the  momentum  equa- 
tions are  given  by  Haney  (1971),  and  only  the  equations  of 
temperature  and  salinity  are  treated  here.   These  are  con- 
veniently written  in  symbolic  form  as 

AS 

-~   =  SALT  ADV.  +  SALT  DIFFUSION  +  CONV.  ADJ.,      (2.9) 

JT  =  TEMP  ADV.  +  TEMP  DIFFUSION  +  CONV.  ADV.,      (2.10) 


where : 


SALT  ADV 


TEMP  ADV 


.'„  -[Vi-lVS)*  ■kUrS)]>  (2.11) 

.  =  -[V*IVT)*  ±l"%  (2.12) 


SALT  DIFFUSION  =   A&V  $  +  *S  A~*J  (2.13) 

TEMP  DIFFUSION  =  ^*^     »  *    "J*L  >  (2.1*0 


and  A$,  fc^,  Aft,  and  *h  are  the  diffusion  coefficients  as 
previously  listed.   Other  numerical  models  have  used  values 
of  rv^lO   cm  sec    and  %**  1  cm  sec"   (Bryan,  1969;  Munk, 
1950).   Because  such  a  simple  form  of  mixing  was  used,  it 
was  reasonable  to  use  similar  values  for  these  coefficients. 
The  coefficients  of  vertical  eddy  diffusion  of  heat  and  salt 
were  fen,k^   =1.0  cm  sec"  .   The  coefficients  of  horizontal 

eddy  diffusion  of  heat  and  salt  were  ^,A§  =  2.5  x  10' 

2    -1 

cm  sec   .   The  friction  terms  in  the  equation  of  motion  were 
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similarly  formulated  in  terms  of  simple  horizontal  and 
vertical  eddy  diffusion  of  momentum  with  coefficients 

^Q  O  ~]  0       1 

¥\=   1.0  x  10   cm  sec   and  V^=   10.0  cm  sec"  .   This  gave  a 

ratio  of  H^  to  K>t  of  ten  which  appeared  to  give  good  results; 

Krivelevich  (1970)  found  a  ratio  of  twenty  to  give  the  best 

results  in  a  forty  layer,  two-dimensional  equatorial  model. 

The  discussion  of  the  "CONVECTIVE  ADJUSTMENT"  is  given 

in  a  later  section,  and  only  the  differencing  of  the  advec- 

tion  and  diffusion  terms  are  treated  here.   Furthermore, 

since  the  forms  of  equations  (2.9)  and  (2.10)  are  identical, 

and  both  variables  S  and  T  are  located  at  the  same  grid 

points,  only  the  finite  differencing  of  the  salinity  equation 

is  given.   However,  because  the  surface  boundary  conditions 

are  different  for  the  vertical  diffusion  terms  of  each,  a 

separate  treatment  of  them  is  necessary.   Although  all 

finite  differencing  was  done  in  curvilinear  coordinates, 

only  the  cartesian  forms  are  shown  for  the  purposes  of 

demonstrating  the  techniques  involved. 

B.   HORIZONTAL  AND  VERTICAL  FINITE  DIFFERENCING 
1.   The  Horizontal  Grid 

A  portion  of  the  horizontal  grid  of  size  21  x  18  is 
shown  in  figure  2.6.   The  temperature  T  and  salinity  S  are 
located  at  the  I  and  ^  grid  points,  whereas  the  velocity 
components  M  and  V  are  located  at  intermediate  points;  ^K 
is  the  E-W  grid  distance,  and  &vj  the  N-S  grid  distance. 
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Figure  2.6   The  horizontal  space  staggered  grid  with 
salinity  S  and  temperature  T  located  at  the  I  and  ^ 
grid  points  and  the  velocity  components  u  and  v  at 
intermediate  points;  fcK  is  the  E-W  grid  distance  and 
tKi    the  N-S  grid  distance. 
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Figure  2.7   Grid  points  near  the  western  wall  used  to 
derive  the  finite  difference  analogues  along  the 
boundary  at  the  general  grid  point  (i,j). 
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2.   Differencing  of  Horizontal  Advection 

Expanding  the  horizontal  advection  portion  of 
equation  (2.11)  gives 


?.(VS)«  ilu^iUS). 


(2.15) 


Considering  first  only  those  grid  points  interior  to  the 
boundaries  and  using  centered  finite  differencing  at  the 
general  grid  point  (i,j)  gives 

7.(V'S)it'  =  (^V^Va.»\"U^l-Vllj  ±   W^lt'a4Vv-CvS)l,>7i.    (2.16) 

Since  neither  u  or  v  is  explicitly  defined  at  the  half 
integer  points,  a  linear  approximation  is  used;  for  example, 
the  term  tu$)i*V^»j   is  approximated  as 


U*S)toYUm  K1^^  4-  U>\iil^ 


2. 


(2.17) 


Using  a  similar  procedure  for  the  remaining  three  terms  in 
equation  (2.16)  gives  the  final  differenced  form  for 
equation  (2.15), 

WVS)i,i-  [mUiw*h+  uu^ii-vj.]  [Slu,j  *  SUJ  /(***) 


"  [>v2.,3-v*.  +^a\,i-vj  [s^-t+S;.j  /(*&)). 


(2.18) 
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The  use  of  a  centered  difference  for  points  along 
the  boundaries  is  not  practical  because  the  velocity  compo- 
nents u  and  v  are  not  specified  exterior  to  the  walls.   The 
finite  difference  analogues  are  derived  by  specifying 
boundary  conditions  and  using  one-sided  approximations  to 
the  derivative  when  necessary.   Because  the  mathematical 
forms  of  the  difference  analogues  are  similar  for  all  four 
walls,  it  is  sufficient  to  demonstrate  the  differencing  for 
only  one  of  them.   A  segment  of  the  western  boundary  region, 
as  pictured  in  figure  2.7,  is  used  to  derive  the  difference 
approximation  at  the  general  point  (l,j). 

Since  the  normal  component  of  the  current  u  must 
vanish  along  the  western  wall,  the  first  term  of  equation 
(2.15)  reduces  to 

A(us)^=   <**>*■* :  ^  (2.i9) 

2. 

Further  reduction   leads   to 

_±(*S)l4-    (  ^VtirVV.  +  >Wi»i«-0(  Sl/i  »  &*.j)   T  (2.20) 

To  derive  the  second  term  of  equation  (2.15),  it  is  assumed 
that  free  slip  applies  for  the  meridional  speed  component  v  , 
that  is, 

Vljj-V^  st   W-lja— *'l   >  and  (2.21) 

VljiW,.  ~   VV^/j^  VL       .  (2.22) 
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In  finite  difference  form  the  term  is 


(2.23) 


Application  of  the  boundary  conditions  given  by  equations 
(2.21)  and  (2.22)  gives 


3.   Differencing  of  Horizontal  Diffusion 

Expansion  of  the  horizontal  diffusion  portion  of 
equation  (2.13)  gives 


A$v**-S  =  As 


(2.25) 


This  is  rewritten  as 


Asv^S  =•  As 


(2.26) 


Evaluation  at  the  general  interior  grid  point  (i,j)  results 
in  the  equation 


M*$L' 


4-  AH~J^V1~7^Kj-^, 


taj 


(2.27) 
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Application  of  a  centered  derivative  to  each  first  order  term 
gives  the  final  result, 


(2.28) 


The  differencing  of  equation  (2.25)  along  the  ocean 
walls  is  handled  similar  to  the  evaluation  of  the  horizontal 
advection  terms  along  the  boundaries;  the  differencing  will 
be  illustrated  for  an  arbitrary  point  (l,j)  along  the 
western  boundary.   Boundary  conditions  along  the  wall  are 
those  of  perfect  insulation, 


4ft  J*"J 


(2.29) 


Writing  equation  (2.25)  in  finite  difference  form  at  point 
(l,j)  gives 


Asv'-Sto 


2. 


(2.30) 


Application  of  the  boundary  conditions  and  further  reduction 
of  the  first  order  terms  gives  the  final  result, 


As^Si^r*  2As 


W1* 


+■  As 


.  (2.3D 


28 


4.   The  Vertical  Grid 

Figure  2.8  shows  the  placement  of  layers  in  the 
vertical.   The  vertical  index  k  and  all  of  the  variables 
(u,  v,  S,  T,  p  )  are  defined  at  the  integer  levels  shown  by 
the  solid  lines,  whereas  Wis  defined  at  the  half  integer 
levels  indicated  by  dotted  lines.   The  (KM  +  1)  unequal 
thicknesses  are  defined  as 


(2.32) 
(2.33) 
(2.3^) 


where  2k<0  while  H>0  is  the  ocean  depth.   This  definition 
satisfies  the  equation 


iH+!^+  i[ 


kt 


&4+y. 


~H. 


(2.35) 


The  thicknesses  centered  at  the  levels  at  which  (u,  v,  S,  T, 
p)  are  located  are  defined  as  the  average  of  the  two  adjacent 
half  integer  levels, that  is, 


^K"  \\    ^*"V2.  *"  ^W-Vl. 


}      K^ljKM. 


(2.36) 


This  definition  has  the  property  that 


(2.37) 
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Figure  2.8   Location  of  the  variables  in  the  vertical  grid 
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5.   Finite  Differencing  of  Vertical  Advection 

Considering  the  intermediate  levels  LtL^.  «ooo  &^fcN\-l> 
the  vertical  advection  portion  of  equation  (2.11)  is  written 
at  the  K^V\  level  as 

*(yrS)<a   (wS^k-Vt  -  (w&Wt  .  (2-38) 

Since  S  is  not  available  at  half  levels,  the  average  of  the 
salinities  at  adjacent  levels  is  used  to  give 

jLtwS}*-   ^k-v».(^i-»SO-vrK»vx(SK»SK40  ,     (239) 

At  the  upper  boundary  the  vertical  velocity  is  zero, 
mathematically 

W-»/2.»  C  (2.40) 


Equation  (2.38)  evaluated  at  K=l  gives 

o 

(2.41) 


±1  &X. 


Further  expansion  gives  the  analogue  at  the  level  nearest 
the  surface  as 

±UrS)L=.   I^^aAJLtJb)  .  (2.42) 
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At  the  ocean  floor  Wkk+^j.  —  0  >  and  at  the  bottom  level 
equation  (2.38)  becomes 

6.   Differencing  of  Vertical  Diffusion 

Considering  the  interior  levels  first,  the  vertical 
diffusion  portion  of  equation  (2.13)  is  written  at  the 
general  k*h  level  as 

c  \ik\        -   4i\   1 

*s^|\  -   ^ftL^jK-Vi    jJ^/jC^L  .  (2.44) 

Further  reduction  of  the  first  order  derivative  terms  gives 

The  lower  boundary  condition  of  perfect  insulation  at 
i.--H   and  index  KM*V2.  is 


J>S  \     _  0  (2.46) 


Using  this  boundary  condition  gives 

-0 


f  ii\       _  iVf      i 


(2.47) 
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and  the  analogue  for  the  lowest  level  becomes 

Up  to  this  point  the  finite  differencing  of  the 
salinity  and  temperature  equations  has  required  only  a 
single  treatment  as  both  equations  have  similar  mathematical 
forms  and  boundary  conditions.   The  difference  analogue 
for  a  term  in  the  heat  equation  is  obtained  from  its  counter- 
part in  the  salinity  equation  by  replacing  S  with  the 
variable  T.   At  the  ocean  surface  this  procedure  does  not 
apply  since  the  boundary  conditions  are  very  different. 

The  upper  boundary  condition  for  the  heat  equation 
is 

kr^l    a  -A.  (2.4q) 

where  Q m  Q  i  +Qi(T^-TlJ  .   Equation  (2.8)  is  a  relatively 
strong  condition  that  ties  Tj^,  the  temperature  at  the 
first  level  below  the  surface,  rather  closely  to  the  speci- 
fied atmospheric  temperature  *X.   Should  some  mechanism 
cause  a  large  difference  between  the  two,  heat  fluxes  are 
created  that  tend  to  restore  the  difference  to  a  small  value 
The  boundary  condition  for  salt  on  the  other  hand  is 


"■SW  SiU'p)- 


(2.50) 
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This  is  a  relatively  weak  condition  of  nearly  constant 
flux  which  tended  to  make  the  model  ocean  respond  very 
slowly  to  the  surface  salt  flux.   Areas  where  the  rate  of 
evaporation  exceeded  that  of  precipitation  tended  to  get 
more  salty  until  advection  and/or  diffusion  became  suffi- 
cient to  create  a  balanced  state.   Similarly,  areas  with 
excess  precipitation  tended  to  become  less  salty  with  the 
approach  toward  steady  state  being  very  slow.   The  reasons 
for  the  slow  salinity  relaxation  time  and  its  effects  upon 
the  dynamics  of  the  model  are  thoroughly  discussed  in 
Section  III. 

The  finite  difference  form  of  the  vertical  diffusion 
term  in  the  salinity  equation  at  level  Kal  is 


HsL3iK~  tiW\ 


f&L 


(2.5D 


Use  of  the  boundary  condition  (equation  2.50)  gives 


(2.52) 


which  reduces   to   the    final    form 


SlIE-?)  fclSi-^) 


Ji^L 


fclrteVa, 


(2.53) 


In  a  similar  manner  the  difference  form  of  the 
vertical  diffusion  term  in  the  heat  equation  at  the  level 


3^ 


k«L  ii 


terA*T\   =        W 


W- 


\ 


lTk)*i.~  "JfcWj 


(2.54) 


&*, 


Using  the  upper  boundary  condition  for  heat  (equation  2.49) 
yields 


(2.55) 


C.   EQUATION  OP  STATE 

1.   Density  as  a  Function  of  Temperature  and  Salinity 

Density  depends  in  a  nonlinear  fashion  upon  tempera- 
ture, salinity,  and  pressure;  however,  the  assumption  of 
incompressibility  reduces  the  dependence  to  two  variables, 
salinity  and  temperature.   If  P0  is  the  density  correspond- 
ing to  arbitrary  reference  values  of  temperature  T&  and 
salinity   So  s    it  follows  from  a  Taylor  series  expansion 
about  To  and  So  that  the  new  value  of  density  at  salinity 
So+  ^  &   an<i  To+fcT  is  given  by 

p(S*+o$,VuT)-=  0($*/Tb)4-£fcU  *i£^T+M*  1  r^+hTiLJpM,  (2.56) 

Data  taken  from  Von  Arx  (1964)  and  reproduced  as 
figure  2.9  shows  that  over  the  range  of  salinities  and 
densities  commonly  found  in  equatorial  oceans,  the  nonlinear 
terms  may  be  dropped  with  little  sacrifice  in  accuracy  and 
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Figure  2.9   Plot  of  sigma-t  vs.  temperature  and  salinity 
after  Von  Arx  (1964) . 
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the  equation  written  as 


^»  fo[  1-AU-To)  V  frli-^jj  (2.57) 

where 

O  •=  the  density  at  temperature  T  and  salinity  S> 

fr0  «  the  density  at  reference  temperature  To  and 
salinity  So^ 

2.   Solutions  for  ot  and  B 

Introducing  the  quantity   Q"^  •=.  (ps,T,c  —  d.)  K  \0 

where  the  subscript  zero  refers  to  sea  level  pressure,  gives 

e^  i-  iSl\  X\tf*  (2.60) 

&a  1^2)        kU>-*  (2.61) 

Po  *S /T«  CONST 


Both  o^  and  £  are  determined  from  figure  2.9  by  graphically 
taking  a  finite  difference  approximation  to  equations  (2 . 60 ) 
and  (2.  6l)  to  give 

^-^L    *."**** -^^V*"**   MS*«"%-*  (2.62) 

pe>  «TT  /sKowrr  P°  (Tc-T&)  "  > 


b»-!-iS:\  iUT*m  Li^-r*)   »\tV  I^K^CLV.1   (2.63) 

The  reference  values  p0,  To,  and  J0  are  1.0276 
gm-cm   ,  278.2  K,  and  3^.72%,  respectively;  this  reference 
state  is  marked  by  a  star  on  figure  2.9.   The  approximations 
for  ^  and  £  are  taken,  however,  at  a  location  on  the  T,S 
diagram  which  is  representative  of  the  temperature  and 
salinity  found  near  the  surface  in  equatorial  oceans  and 
the  dependence  of  density  on  these  variables  is  nearly 
linear. 

D.   INSTANTANEOUS  CONVECTIVE  ADJUSTMENT  OF  DENSITY 

Whenever  a  heavier  water  mass  lies  above  one  of  lesser 
density,  an  unstable  situation  exists  and  natural  buoyancy 
forces  cause  vertical  mixing  and  a  return  to  a  neutral 
state.   Since  from  the  equation  of  state,  density  is  a 
function  of  temperature  and  salinity,  any  mechanism  that 
causes  temperature  to  increase  or  salinity  to  decrease  with 
depth  is  a  potential  source  of  instability.   In  this  model 
the  only  heat  flux  into  the  ocean  occurs  at  the  surface.   In 
a  steady  state,  the  surface  averaged  heat  flux  must  be  zero. 
Thus  there  are  some  regions  where  the  net  heat  flux  at  the 
surface  is  downward  and  some  regions  where  it  is  upward, 
and  it  is  in  these  regions  of  upward  heat  flux  that  thermally- 
induced  instability  is  found.   Salinity-induced  instability 
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occurs  in  regions  of  the  model  where  evaporation  exceeds 
precipitation  causing  a  layer  of  heavy,  salty  water  to 
exist  near  the  surface.   The  vertical  juxtaposition  due  to 
advection  of  water  masses  of  equal  temperature  but  differ- 
ent salt  content  is  also  a  potential  source  of  salinity- 
induced  instability. 

Instability  in  the  real  ocean  lasts  only  a  short  period 
of  time.   Buoyancy  forces  produce  large  vertical  velocities 
which  mix  the  unstable  layers  very  rapidly  and  efficiently, 
and  quickly  establish  a  uniform  vertical  distribution  of 
mixed  quantities.   The  following  is  a  reasonable  method  of 
modeling  this  phenomena. 

Figure  2.10  shows  a  vertical  section  of  the  ocean  where 
the  density  increases  with  depth  and  represents  the  stable 
case  where  no  convective  instability  will  exist.   Using 
previously  introduced  notation,  several  levels  K-l,  K,  and 
K*l.  and  the  vertical  thickness  ^"^fe  are  shown.   The  inter- 
mediate levels  K-3/2  ,  K.-V2.  ,   K+y^  and  K4-3/^  are  shown  as 
dashed  lines  and  thicknesses  corresponding  to  fc"  «£,  and 
K*1/^  are  denoted  as  ^Ak-Vi.  and  AtM-'/j,  respectively . 
Although  in  the  real  ocean  the  profile  can  be  nonlinear,  a 
linear  profile  as  denoted  by  the  solid  lines  in  figure  2.11 
is  assumed  in  the  model;  this  first  order  approximation  is 
consistent  with  others  used  elsewhere  and  is  sufficient  to 
guarantee  conservation  of  the  quadratic  quantity  in  the 
vertical  advection  terms.   Additionally,  it  is  assumed  that 
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Figure  2.10   Stable  Density  Profile 


K-l 


K 


K-H 


Figure  2.11   Model  Equivalent  of  the  Stable  Density  Profile 


40 


the  density  which  is  stored  at  each  level  represents  the 
vertically  integrated  average  of  the  entire  fc^+h  layer. 

If  Pk-1>Pk   f°r  anv  two  consecutive  levels,  then 
convective  instability  on  a  scale  which  is  too  small  to 
be  resolved  by  the  grid  will  occur.   Since  from  equation 
(2.57)  density  is  a  function  of  both  salinity  and  tempera- 
ture, either  one  or  both  may  be  the  cause  of  the  instability 

Figure  2.12  shows  the  case  where  density  instability 
is  due  to  temperature;  salinity  is  assumed  to  be  constant 
with  depth.  "Tk-i.  represents  the  temperature  for  the  upper 
layer  of  thickness  ^K-t  and  T^  the  temperature  for  the 
lower  layer  of  thickness  *&*,   .   To  insure  stability,  new 
temperatures  T^-i  and  T*  are  determined  such  that  a 
slightly  stable  temperature  lapse  rate  results.   Mathemati- 
cally the  following  must  hold: 

where  &t  is  a  small  positive  constant  whose  value  is  the 
new  lapse  rate  shown  as  the  short  dotted  line  in  figure 
2.12.   To  insure  conservation  of  heat,  the  thickness 
weighted  average  of  the  temperatures  before  and  after  mixing 
must  be  the  same,  that  is, 


"Tk-L  ^k-1  +  T*  fh K  *  "T*£  AK.L  +  Tk* A*K  .       (2.65) 
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Figure  2.12   Diagram  of  an  unstable  state  due  to  the 
temperature  profile  showing  the  necessary  temperature 
adjustments  to  achieve  stability  between  levels  %■*.   and 
2k.-1-   The  solid  line  is  before  adjustment,  the  dashed 
is  after. 
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Equations  (2.6*0  and  (2.65)  comprise  a  set  of  two  simultane- 
ous equations  in  two  unknowns  and  when  solved  give 

T*«   pQjg-K  +  (  Tk-i  -  er&%K-Vi)  ^K-t  ^         (2.66) 

TK_*  a  Tk  4-  €T&%K-ya..  (2.67) 

Equations  (2.66)  and  (2.67)  define  the  new  temperature  values 
If  the  temperature  is  constant  but  a  decrease  in  salinity 
with  depth  causes  instability,  convective  adjustment  must 
take  place  such  that  the  salinity  values  S^  and  S<-jl  are 
replaced  by  new  ones  S^  and  $*«!  such  that  a  stable, 
slightly  positive  vertical  salinity  gradient,  &$ ,  is 
established.   Assuming  conservation  of  salinity,  and  using 
a  reasoning  similar  to  that  used  for  temperature,  gives  the 
following  results  analogous  to  (2.66)  and  (2.67): 

$K*~    Sk&^k  *  (  Sk-i  4-  es&^k-vQ^K-i  i  (2.68) 


S*i*m   SK  -  fc^fcte-Va..  (2.69) 


For  an  incompressible  ocean,  the  value  for  £«r  is  zero, 
but  it  is  reasonable  to  define  it  to  be  some  small  positive 
value  to  avoid  possible  instability  due  to  round  off  error 
in  the  computation;  ^t  is  defined  in  terms  of  To,  the 
surface  reference  temperature,  and  H,  the  depth  of  the  model 
ocean,  as 
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H 


(2.70) 


Having  fixed  £?,  the  appropriate  value  for  fc$  and  €p 

can  be  derived  from  the  equation  of  state  and  the  definitions 
of  €$  (2.64)  and  6T(2.70).   If  it  is  assumed  that  the 
salinity  and  temperature  have  equal  weights  in  determining 
the  density  lapse  rate,  then  it  can  be  shown  from  the  equa- 
tion of  state  that 


fc$=  sL  eTj  (2.71) 


$?  m  2<X€t.  (2.72) 

A  very  brief  explanation  of  how  the  convective  adjustment 
is  programmed  is  as  follows: 

(1)  The  density  lapse  rate,  ((,  ^>K-i  -  ^  /  &£k.-vO  ,  is 
computed  for  a  layer  and  compared  with  the  stable  value  (f  . 

(2)  If  the  density  lapse  rate  is  more  stable  than  Cv>,  no 
change  is  made  in  temperature  or  salinity;  if  it  is  less 
stable  than  €j>,  the  temperature  and  salinity  are  both 
adjusted  using  equations  (2.66),  (2.67),  (2.68),  and  (2.69). 

(3)  All  layers  are  processed  in  order  from  the  first  below 
the  surface  to  the  last  above  the  lower  boundary.   In  the 
convective  adjustment  of  a  pair  of  levels,  it  is  possible 
to  render  an  adjacent  pair  of  levels  unstable;  therefore, 
the  adjustment  proceeds  in  repeated  scans  until  no 


Hk 


adjustment  of  temperature  or  salinity  is  performed  during  a 
scan  (see  figure  2.12).   The  convective  adjustment  is  then 
terminated. 
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III.   RESULTS 

A.   THE  TIME  DEPENDENT  BEHAVIOR  OP  THE  MODEL 
1.   Salinity  Relaxation  Time 

The  upper  boundary  condition  for  the  salinity  equa- 
tion was  the  flux  of  salt  into  an  initially  isohaline  ocean. 
However,  some  insight  into  the  time  dependent  behavior  of 
the  model  can  be  gained  from  the  analysis  of  a  more  idealized 
situation,  one  where  the  upper  layer  of  an  initially  salt- 
free  ocean  of  infinite  depth  is  raised  to  a  constant  value 
S.   Neglecting  the  effects  of  advection,  horizontal  diffu- 
sion, and  convective  adjustment  the  salinity  equation  (2.9) 
and  its  initial  and  boundary  conditions  become 


iL  =  KS^S  (3.D 

S(t,o)=-  Constat*  $*  ^or  t>0,  (3.2) 


S(0,*)«  O.  (3.3) 

To  solve  equation  (3.1)  two  new  variables  T  and  g"  are 
defined  as 

X-  kst,  (3.4) 


tr-  -?—  =    *       .  (3.5) 
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Transforming  (3.1)  into  T  and  <T  space  gives 

il«   '  **  .  (3.6) 

Assuming  that  &  depends  on  time  t  only  through  (j*  implies 
that 


(3.7) 


or  that 


4i  =.  z£.  |**\ 


(3.8) 


Substitution  of  this  result  into  equation  (3.6)  gives 


iii-  =.  -l^iL.  (3.9) 

Letting  Q.-  ^  yields 

^>  =  -2L(T<V)  (3.10) 


whose    solution   is   well   known   as 


%m   %oCTrX.  (3.1D 


Further  integration  of  (3-9)  gives 

o       » 
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(3.12) 


where 

S*~  &(f»°)»  S(,sa»0)i  t>0.  (3.13) 

Finally  the  solution  becomes 

Stir)  -  S*  v  )vc  4*  (3.1^) 

where  S  is  the  value  of  S  at  ^«0. 

The  value  of  <ko  is  determined  from  equation  (3.3)  as  follows 

$(&,■&)  -  SC^)-0  for  fc?0,  (3.15) 

from  equation  (3-1*0 

$i»)  =   $V  ^  $  t^^o,  (3.16) 

since   IcT1  A?  *  ^/Z>  (3.17) 

the  final  solution  for  Q*  becomes 

V  J1F=-  S*  ^3.18) 

and  the  solution  of  equation  (3.1)  including  the  boundary 
condition  becomes 

$(<r)=-  %*l\-  £^(<r))  (3.19) 
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where 


ERKri*   h.  \  eT1  *«?.  (3.20) 

Figure  3.1  is  a  plot  of  equation  (3.19)  and  was 
used  to  determine  an  estimate  for  the  salinity  relaxation 
time.   From  this  graph  it  was  clear  that  fl*»0  gives  Sa  Sj 
the  true  steady  state  solution;  however,  since  both  the 
depth,  £  ,  and  the  vertical  diffusion  coefficient  for  salt, 
fc$j  were  nonzero,  a  true  steady  state  would  exist  only  as 
time  approached  infinity.   To  avoid  ambiguity,  it  was 
convenient  to  define  a  quasi-steady  state  to  exist  when 
S^O.^S   which  corresponded  to  a  j  value  of  0.1.   Solution 
of  equation  (3.5)  for  *t  gives  the  time  required  for  the 
level  £  to  be  affected, 

*«  tf*  r  (3.21) 

Since  *t  varied  directly  as  the  square  of  the  depth, 
the  salinity  relaxation  time  was  expected  to  be  very  long 
for  the  lower  layers.   For  ^g  a  0.1  cm  sec"   and  f"  =  0.1 
the  first  layer  at  five  meters  computed  to  be  steady  state 
in  0.2  years,  the  third  in  20  years,  the  fifth  in  500  years, 
and  the  seventh  in  32,000  years.   As  the  model  was  inte- 
grated over  the  relatively  short  period  of  30  years,  theo- 
retically only  the  first  90  meters  were  expected  to  be  in 
equilibrium.   Furthermore,  since  the  upper  boundary  condition 
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Figure  3.1  A  plot  of  S*  vs.  fal/W);  the  dotted 
line  shows  a  quasi-steady  state  value  of   0.<\S*  . 
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was  the  specification  of  the  flux  of  salt  and  not  the 
idealized  condition  treated  above,  the  relaxation  time  was 
expected  to  be  even  larger  for  each  level  than  indicated 
by  (3.21). 

Figure  3.2  is  the  salinity  anomaly  vs.  time  plot 
at  various  levels  for  a  30-year  period  at  the  grid  point 
X=.lOj  T*  B;   (see  figure  3.20).   This  location  was  chosen 
for  studying  the  dependence  of  the  salinity  relaxation  time 
upon  the  vertical  diffusion  because  the  prescribed  upward 
flux  of  salt  at  the  surface  was  a  maximum  here.   Addition- 
ally, because  both  the  horizontal  and  vertical  currents 
were  small  with  little  shear,  the  effects  of  advection  were 
expected  to  be  negligible.   The  apparent  breaks  at  approxi- 
mately 10  years  in  the  otherwise  smooth  curves  were  caused 
by  a  change  in  h|^,  the  vertical  diffusion  coefficient  for 
momentum,  and  did  not  affect  the  general  nature  of  the 
adjustment.   As  predicted  from  the  theoretical  analysis, 
levels  above  100  meters  appeared  to  be  approaching  steady 
state  and  showed  little  change  during  the  last  10  years  of 
integration;  however,  levels  below  250  meters  showed  a 
constant  rate  of  increase  of  approximately  0.002*/oo  per 
century  at  the  end  of  the  integration.   Although  the  salinity 
at  the  lower  levels  was  not  in  steady  state,  it  was  felt  that 
this  was  not  a  serious  problem  because  the  salinity  gradient 
did  not  have  a  significant  effect  upon  the  model  dynamics, 
which  were  primarily  controlled  by  the  surface  wind  stress 
and  the  near  surface  eddy  momentum  transfer  processes. 
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Figure  3.2   A  plot  of  salinity  vs.  time  at  the  grid  point 
with  maximum  upward  surface  flux  for  levels  5,  20,  100, 
250,  and  700  meters.   The  ordinate  values  have  been  shifted 
to  place  all  curves  on  one  graph  and  only  the  slopes  are 
significant . 
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2.   Temperature  Relaxation  Time 

As  discussed  in  Section  II,  the  upper  boundary 
condition  for  temperature  was  quite  different  than  the  one 
for  salinity.   The  flux  of  salt  at  the  ocean  surface  was 
specified  as  a  function  of  the  existing  salinity  and  the 
difference  between  the  fixed  rates  of  evaporation  and  pre- 
cipitation (equation  2.50).   If  there  were  no  horizontal 
transport  of  salt,  a  steady  state  could  not  be  reached. 
Because  the  prescribed  flux  is  nearly  independent  of  salin- 
ity, those  regions  with  precipitation  excesses  would  tend 
to  become  less  salty,  whereas  those  with  net  deficits  would 
tend  to  become  more  salty.   Temperature,  on  the  other  hand, 
was  different  because  the  ocean  surface  temperature  was 
tied  to  the  apparent  atmospheric  equilibrium  temperature 
Tfc*  es.   TA  +  &1/&2.  (equation  2.M9).   If  the  temperature  of 
the  near  surface  layer  T^  became  much  different  than   «X  j 
then  a  large  heat  flux  resulted  which  drove  Ti  towards  *h 
and  resulted  in  a  reduction  of  the  heat  flux.   If  there  were 
no  horizontal  transport  of 'heat,  a  steady  state  could  even- 
tually exist  in  which  the  temperature  of  the  ocean  was  every- 
where  equal  to  l/^ ,    For  these  reasons,  it  was  apparent 
that  the  theoretical  analysis  just  performed  applied  more 
closely  to  temperature  than  to  salinity.   Accordingly,  it 
was  expected  that  the  time  required  to  reach  temperature 
equilibrium  would  be  less  than  that  needed  for  salinity  to 
reach  steady  state.   To  decrease  the  temperature  relaxation 
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o 
time,  an  initial  temperature  lapse  rate  of  -1  /100  meters 

was  chosen,  a  value  that  was  approximately  equal  to  the 

steady  state  lapse  rate  in  the  equatorial  domain  of  Haney's 

(1971)  model. 

Figure  3.3  shows  the  temperature  vs.  time  plots  for 

several  levels  at  the  same  point  as  in  figure  3.2.   As 

expected,  the  near  surface  levels  reached  steady  state 

comparatively  quickly,  whereas  the  deeper  ones  did  not  reach 

equilibrium  over  the  30-year  integration  period.   The  upper 

four  levels  achieved  equilibrium  after  only  10  years  and 

then  showed  little  change  for  the  next  20  years,  whereas 

the  corresponding  salinity  values,  although  approximately 

steady  state,  still  continued  to  show  some  rate  of  change 

over  the  corresponding  period.   At  the  end  of  30  years,  the 

bottom  ocean  levels  appeared  to  be  warming  at  a  constant 

o 
rate  of  0.25  C  per  century.   This  same  trend  was  noted  in 

Bryan's  model  (1968)  although  his  time  integration  was 

carried  out  over  a  much  longer  time  span.   As  in  Bryan's 

model,  the  gradual  heating  was  probably  due  to  the  selection 

of  too  cold  of  initial  temperatures  for  the  lower  layers. 

Again  as  for  salinity,  the  failure  to  reach  steady  state 

in  the  bottom  layers  was  not  felt  to  be  especially  damaging 

to  the  dynamical  results  obtained.   The  circulation  was 

determined  to  be  primarily  wind-driven,  and  it  is  the 

surface  layers  which  are  of  major  interest  in  an  ocean 

forecasting  model. 
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Figure  3.3   Temperature  vs.  time  plots  for  levels  5,  20,  100, 
250,  and  700  meters  over  a  30-year  period.   The  ordinate 
values  have  been  transformed  various  amounts  to  place  all 
curves  on  one  plot  and  only  the  slopes  are  significant. 
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3.   Dynamic  Current  Relaxation  Time 

Munk  (1950),  Charney  (1959),  Veronis  (1959),  and 
other  researchers  have  proved  that  the  wind  plays  an  impor- 
tant role  in  determining  the  large  scale  features  of  the 
ocean  circulation.   It  is  generally  known  that  the  negative 
zonal  wind  stress  in  the  near  equatorial  region  causes  a 
strong  surface  current  towards  the  west,  and  that  Ekman 
divergence  is  responsible  for  a  meridional  component  of 
circulation  at  the  surface  towards  each  pole.   Additionally, 
it  has  been  established  that  the  ocean  response  time  to 
the  surface  stresses  is  very  short,  approximately  100  to 
200  days. 

No  formal  mathematical  analysis  of  the  current 
response  time  is  attempted  here,  and  only  the  actual  results 
of  the  time  integration  are  presented  to  verify  that  the 
dynamic  response  time  is  small.   Figures  3.  4,  3-5,  and  3-6 
are  plots  of  the  zonal  current  taken  at  the  equator  midway 
between  the  eastern  and  western  boundaries.   The  computer 
analyzed  data  was  taken  from  the  time  integration  every  25 
days.   Figure  3-^  is  representative  of  the  second  level 
below  the  surface  where  a  strong  countercurrent  toward  the 
east  of  180  cm/sec  had  developed  during  the  first  10  years. 
Figure  3.5  shows  the  flow  at  100  meters  where  a  sizeable 
countercurrent  continued  to  exist  after  adjustment  of  fc^. 
Figure  3.6  is  the  zonal  current  vs.  time  at  700  meters  and 
represents  typical  deep  ocean  water.   The  response  time  of 
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Figure  3«^   A  plot  of  zonal  current  vs.  time  at  the  grid 
point  I—  10,  T«*  S   located  on  the  equator  midway  between 
the  eastern  and  western  boundaries.   The  depth  is  20 
meters . 
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Figure  3-5   Zonal  current  u  (cm/sec)  vs.  time  for 
grid  location  J>  \0 s   l*°l     at  100  meters. 
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Figure    3.6      Zonal   current   u  (cm/sec)    vs.    time   at   grid 
point     3>10,    3=1    at    700   meters. 
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about  100  days  was  extremely  fast  compared  to  those  of 
salinity  and  temperature  and  appeared  to  be  Independent  of 
depth. 

The  rapidity  with  which  the  wind-driven  circulation 
adjusts  was  dramatically  demonstrated  after  h^  was  altered 
to  correct  two  unrealistic  features  of  the  circulation  that 
had  appeared  during  the  first  10  years.   An  unreasonably 
strong  equatorial  surface  current  had  developed,  and  its 
associated  countercurrent  was  located  only  20  meters  below 
sea  level.   Within  100  days  after  increasing  kw,  the  currents 
at  most  levels  had  adjusted  to  new  values  radically  differ- 
ent from  the  previous  values.   At  20  meters  the  flow  shifted 
from  60  cm/sec  from  the  west  to  kf   cm/sec  from  the  east. 
Additionally,  the  core  of  the  countercurrent  now  appeared 
at  100  meters,  a  depth  more  representative  of  its  position 
in  the  real  ocean. 

B.   EQUATORIAL  UNDERCURRENT 

One  of  the  unique  features  of  the  Pacific  Ocean  is  the 
existence  of  a  thin,  swift  countercurrent  which  is  symmetric 
about  the  equator  and  moves  eastward  below  the  westward 
flowing  South  Equatorial  Current.   This  phenomena,  called 
the  Cromwell  Current,  was  first  studied  in  detail  by  Knauss 
(I960)  who  determined  the  current  to  be  about  300  km  wide, 
approximately  200  ra  thick,  and  with  velocities  as  high  as 
150  cm/sec.   The  core,  a  jet-like  structure  with  stronger 
shear  above  its  center  than  below  it,  was  found  to  extend 
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over  at  least  3500  n  mi  at  a  depth  that  varied  from  100 
meters  at  140  W,  to  about  75  meters  near  the  Galapagos 
Islands.   In  a  later  study  Knauss  (1966)  found  a  much  weaker 
flow,  and  concluded  that  the  countercurrent  probably  exhib- 
ited considerable  year-to-year  variability.   He  theorized 
that  as  a  first  approximation,  the  countercurrent  could  be 
considered  as  a  separate  flow  imbedded  in  the  surface  South 
Equatorial  Current.   A  stronger  than  average  surface  current 
was  thought  to  be  accompanied  by  a  weaker  undercurrent, 
whereas  a  weaker  surface  current  was  believed  to  be 
associated  with  a  stronger  than  average  undercurrent. 

Several  theoretical  analyses  of  the  Cromwell  Current 
have  been  performed.   Charney  (1959)  showed  that  for  a 
homogeneous  ocean  driven  by  a  uniform  wind,  the  zonal  current 
profile  was  parabolic  with  the  maximum  velocity  dependent 
upon  the  coefficient  of  vertical  eddy  viscosity  v>  .   For 
large  values  of  u  the  easterly  surface  current  extended  to 
the  hypothetical  level  of  no  motion.   As  V  was  allowed  to 
decrease,  the  countercurrent  began  to  develop;  however,  at 
very  small  values  of  U  ,  the  countercurrent  became  so  well 
developed  that  the  flow  at  the  surface  had  reversed  and 
become  westerly.   Veronis  (1959)  offered  a  simplified 
explanation  of  the  undercurrent  using  the  idealized  basin 
shown  as  figure  3.7.   He  reasoned  that  a  surface  wind  blowing 
from  the  east  caused  an  accumulation  of  water  at  the  western 
boundary,  which  caused  the  pressure  gradient  to  drive  the 
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Figure  3.7   (After  Veronis)   A  cellular  circulation 
develops  when  a  nonrotating  basin  is  subjected  to 
a  surface  wind  stress;  the  increased  pressure  at  -  the 
western  boundary  causes  an  eastward  countercurrent 
to  exist  below  the  westward  surface  current. 
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subsurface  water  towards  the  east.   Arthur  (1959),  by 
neglecting  all  terms  in  the  x-component  equation  at  the 
equator  except  the  pressure  gradient  and  frictional  stress 
ones,  showed  that  a  representative  zonal  profile  containing 
a  countercurrent  could  be  calculated  from  values  of  surface 
current  velocity,  wind  stress,  and  pressure  gradient  for  a 
constant  U  .   Additionally,  he  concluded  that  although  the 
solution  was  stable,  the  nonlinear  terms  win  and  W^&were 
probably  not  negligible.   Unfortunately,  inclusion  of  the 
nonlinear  terms  would  have  made  the  analytical  analysis 
difficult. 

In  an  equatorial  prediction  model  with  good  vertical 
resolution  near  the  surface,  the  equatorial  countercurrent 
should  be  readily  resolved.   Figure  3.8  is  an  E-W  cross 
section  of  zonal  current  taken  along  the  model  equator  at 
the  end  of  the  30  years  of  time  integration.   A  counter- 
current  of  about  20  cm/sec  had  developed  at  a  mean  depth  of 
about  100  meters.   This  current  appeared  to  be  slightly  more 
intense  and  at  a  greater  depth  near  the  western  boundary 
than  near  the  eastern  boundary,  and  exhibited  considerably 
more  vertical  current  shear  above  the  core  axis  than  below 
it.   Additionally,  at  approximately  700  meters  the  flow  was 
toward  the  west,  as  it  was  at  the  surface.   For  comparison, 
the  corresponding  equatorial  cross  section  for  the  Pacific 
after  Knauss  (I960)  is  shown  as  figure  3.9.   Although  the 
observed  data  exhibited  considerable  station  to  station 
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variability,  all  of  the  gross  features  are  easily  identified 
in  the  model.   However,  significant  differences  between  the 
observational  and  model  results  were  found  to  exist.   The 
velocity  of  the  Cromwell  Current  averaged  about  140  cm/sec 
compared  to  the  model  mean  value  of  20  cm/sec.   Also,  the 
mean  upward  slope  of  the  current  core  toward  the  east  was 
somewhat  greater  in  the  Pacific  cross  section;  this  may  be 
due  to  insufficient  vertical  resolution  in  the  model. 

As  a  further  comparison,  the  vertical  profile  of  the 
east-west  averaged  current  in  the  model  was  plotted,  (figure 
3.10),  with  a' typical  Pacific  profile  determined  by  Knauss 
(1966).   Although  the  magnitudes  were  quite  different,  the 
vertical  current  shear  above  the  core  of  the  undercurrent 
was  quite  similar,  and  both  were  approximately  parabolic 
in  agreement  with  Charney's  theoretical  results.   Below  the 
point  of  maximum  velocity,  the  real  data  retained  its  para- 
bolic shape,  whereas  the  mathematical  form  of  the  model 
profile  is  not  easily  determined.   One  possible  explanation 
for  the  difference  in  the  profiles  is  that  the  mean  surface 
winds  during  the  observation  period  were  below  average. 
Then  according  to  Knauss' s  concept  of  the  undercurrent  as 
a  flow  embedded  in  the  surface  current ,  the  weaker-than- 
normal  surface  winds  would  have  caused  a  below  average 
surface  flow  and  an  above  average  countercurrent .   Although 
this  simple  argument  has  some  merit,  a  more  likely  explana- 
tion lies  in  the  way  vertical  diffusion  of  momentum  was 
handled  in  the  model.   The  assumption  of  a  constant  kinematic 
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Figure  3.10   Plot  comparing  the  average  model  zonal  current 
shown  by  the  dashed  line  with  the  smoothed  zonal  current  at 
118°W  in  the  eastern  tropical  Pacific  as  determined  by 
Knauss  (1966).   The  model  time  is  30  years. 
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eddy  viscosity,  Kmj  implied  that  the  eddies  that  transfer 
momentum  vertically  are  independent  of  depth.   In  the  near 
surface  layers  of  the  ocean,  *k  may  increase  toward  the 
surface  due  to  increased  wind  mixing,  but  the  exact  func- 
tional dependence  of  ^m  on  £  is  not  clear.   Also,  Km  should 
depend  upon  the  vertical  stability  of  the  water  column  as 
measured  perhaps  by  a  large  scale  Richardson  number,. but 
again  the  best  formulation  of  the  dependence  is  not  known. 

The  jet-like  structure  of  the  model  undercurrent  was 
well  displayed  in  a  N-S  cross  section  of  zonal  current 
taken  approximately  midway  between  the  west  and  east  bound- 
aries  (figure  3.11).   The  strong  wind-generated  surface 
current,  its  associated  undercurrent,  and  the  pronounced 
vertical  current  shear  were  all  readily  apparent.   These 
significant  features  were  also  apparent  in  figure  3.12, 
a  typical  cross  section  taken  by  Knauss  on  his  1958  expedi- 
tion, during  which  the  current  was  the  stronger  of  the  two 
expeditions.   Also  shown  on  figure  3.12  is  the  thermal 
structure  of  the  undercurrent.   Although  the  isotherms 
showed  considerable  small  scale  spatial  variations,  the  large 
scale  pattern  revealed  that  the  isopleths  were  concave 
downward  above  the  current  core,  whereas  below  it  the  con- 
cavity was  reversed.   In  an  attempt  to  identify  a  similar 
thermal  structure  in  the  model,  the  model  countercurrent 
was  superimposed  on  the  N-S  cross  section  of  temperature 
(figure  3.13).   A  cold  region  above  the  undercurrent  axis 
was  evidenced  by  isotherms  that  were  concave  downward. 
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However,  this  same  pattern  persisted  for  200  meters  below 
the  axis,  where  only  a  slight  reversal  in  curvature  was 
noted.   The  failure  of  the  model  to  adequately  resolve  this 
thermal  structure  is  most  likely  due  to  the  simple  form  of 
the  kinematic  eddy  viscosity.   The  thermal  structure  in 
the  vicinity  of  the  equator  in  figure  3.13  is  adequately 
explained  by  vertical  advection  by  a  strong  upward  vertical 
velocity  above  about  300  m,  and  weaker  downward  motion  below 
300  m  (figure  3.23).   The  equatorial  isotherm  spreading  in 
the  real  ocean  (figure  3.12)  is  probably  caused  by  enhanced 
vertical  mixing  in  the  strong  shear  regions  above  and  below 
the  jet  core.   An  important  goal  of  any  future  research  on 
this  model  should  be  to  determine  a  more  adequate  vertical 
eddy  diffusion  coefficient,  Ki-\  ;  one  that  is  a  function  of 
vertical  stability,  vertical  current  shear,  and  depth. 

C.   CURRENT  CIRCULATIONS 

The  model  was  driven  primarily  by  the  prescribed  surface 
wind  stress  T^  shown  in  figure  2.3.   The  existence  of  the 
single  minimum  of  Tfc  displaced  slightly  north  of  the  equator 
is  important,  for  it  causes  the  curl  of  the  wind  stress  to 
be  greater  in  the  northern  domain.   As  Munk  (1959)  showed, 
this  asymmetry  would  tend  to  cause  a  stronger  circulation 
in  the  Northern  Hemisphere,  and  the  non-zero  wind  stress 
curl  on  the  equator  would  cause  a  circulation  gyre  to  exist 
there.   To  more  fully  discuss  this  last  point  the  general 
circulation  of  the  model  as  a  whole  will  be  described. 
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Figures  3.1*1  and  3.15  are  the  zonal  and  meridional 
current  components  at  five  meters  and  illustrate  the  general 
flow  patterns.   Near  the  equator  a  strong  zonal  current 
flow  towards  the  west  had  been  created.   Because  the  coriolis 
force  is  exactly  zero  on  the  equator,  there  was  little  tend- 
ency  for  this  current  to  flow  laterally  away  from  the 
equator  and  the  meridional  current  component  was  essentially 
zero.   Several  degrees  on  each  side  of  the  equator  a  geo- 
strophic  plus  friction  balance  tended  to  take  place  and 
Ekman  divergence  was  evident.   The  zonal  component  decreased 
in  magnitude  and  the  meridional  current  increased  as  the 
surface  mass  transport  tended  to  become  essentially  poleward 
in  both  northern  and  southern  hemispheres.   In  response  to 
the  divergence  of  surface  water,  vertical  upwelling  of  cool 
bottom  water  took  place  in  the  equatorial  region  to  maintain 
mass  continuity.   The  weak  southward  directed  flow  on  the 
equator  in  the  eastern  half  of  the  basin,  and  alternating 
north-south  flow  in  the  western  half,  is  a  part  of  the 
horizontal  gyre  driven  by  the  curl  of  the  wind  stresses 
mentioned  above. 

Figure  3.16,  a  horizontal  cross  section  of  vertical 
velocity,  w",  showed  the  existence  of  a  strong  band  of  upwell- 
ing centered  on  the  equator.   This  figure  also  showed  the 
effects  of  the  boundaries  on  the  flow  pattern.   Along  the 
eastern  boundary,  strong  upwelling  was  necessary  to  replace 
the  surface  water  carried  west  by  the  strong  surface  current; 
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near  the  western  boundary  the  large  accumulation  of  surface 
water  due  to  the  surface  current  caused  large  downwelling. 
Near  both  northern  and  southern  boundaries,  strong  down- 
welling  was  necessary  to  compensate  for  the  poleward  trans- 
port of  surface  water  in  both  hemispheres.   The  nonuniform 
structure  of  the  horizontal  cross  section  of  meridional 
current  (figure  3.15)  is  caused  by  the  presence  of  Rossby 
waves  in  the  vertical  mean  part  of  the  current  that  have 
propagated  towards  the  west,  impinged  upon  the  western 
boundary,  and  become  standing  waves  (Munk,  1950). 

D.   THERMAL  STRUCTURE 

Since  the  model  was  driven  primarily  by  the  surface 
wind  stress,  it  was  anticipated  that  the  patterns  of  current 
flow  would  have  a  large  effect  upon  the  thermal  structure. 
Figure  3.17,  a  horizontal  cross  section  of  temperature  at 
five  meters,  showed  this  to  be  the  case.   A  relatively  cold 
tongue  caused  by  upwelling  was  formed  at  the  equator.   This 
minimum  was  most  pronounced  at  the  eastern  boundary  where 
cold  water  was  advected  to  the  surface  by  the  strongest 
vertical  velocities.   As  this  cooler  water  was  transported 
to  the  west,  it  gained  heat  and  the  thermal  gradient  weak- 
ened.  Two  warm  tongues  were  formed  at  6  N  and  9  S  due  to 
surface  heating  and  weak  sinking  of  the  cool  equatorial 
water  after  it  had  diverged  from  the  equator.   Further  pole- 
ward movement  of  this  equatorial  water  caused  cooling  as 
the  amount  of  surface  heating  began  to  decrease. 
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A  realistic  modeling  of  the  vertical  thermal  structure 
is  essential  for  any  ocean  model  that  is  to  be  coupled  to 
an  atmospheric  model  for  long  range  prediction  purposes. 
Figure  3«l8  is  an  E-W  vertical  cross  section  of  temperature 
along  the  equator  and  showed  a  thermal  gradient  that  was 
gravitationally  stable.   Except  for  the  layers  above  100 
meters,  which  exhibited  a  nearly  neutral  temperature  gradi- 
ent, all  layers  showed  a  large  stable  lapse  rate.   The  ini- 
tial  temperature  gradient  of  -1  C  per  100  meters  had  not 
changed  appreciably  over  the  30  year  time  integration  for  levels 
below  200  meters.   What  appeared  to  be  the  top  of  the  perma- 
nent thermocline  existed  at  depths  that  varied  from  almost 
zero  near  the  eastern  boundary  to  about  100  meters  near  the 
western  wall.   The  near  isothermal  temperatures  normally 
found  in  real  oceans  at  great  depths  had  not  developed 

appreciably.   However,  towards  the  end  of  the  integration 

t 
the  two  lower  levels  were  warming  at  the  rate  of  0.25  C  per 

decade,  and  it  was  felt  that  over  a  very  long  time  period 
the  desired  result  could  have  been  realized. 

To  further  examine  the  thermal  structure,  vertical  pro- 
files of  temperature  vs.  depth  at  the  three  selected  stations 
marked  on  figure  3.18  are  shown  as  figure  3-19.   The  first 
station  at  10°  west  of  the  eastern  boundary  showed  no  iso- 
thermal layer;  equatorial  upwelling  had  tended  to  force  the 
characteristics  of  the  deep  water  on  the  surface  layers. 
The  next  station  located  at  the  midpoint  between  the  bound- 
aries showed  the  development  of  a  near  isothermal  layer  due 
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to  the  shallower  vertical  velocity  pattern  toward  the  west. 
Station  number  three  located  10°  east  of  the  western  boundary 
showed  the  deepest  isothermal  layer.   All  three  profiles  had 
a  significant  gradient  between  layers  one  and  two,  a  distance 
of  only  fifteen  meters.   In  the  real  ocean  a  weak  wind  would 
have  been  sufficient  to  mix  this  layer  and  establish  iso- 
thermal conditions.   One  of  the  weak  points  of  the  present 
model  is  felt  to  be  the  lack  of  an  effective  modeling  of 
surface  wind  mixing  and  internal  vertical  mixing  in  the 
stable  case. 

E.   SALINITY  STRUCTURE 

The  salinity  distribution  was  largely  determined  by  the 
upper  boundary  condition  for  salt  which  depended  upon  the 
net  rates  of  evaporation  and  precipitation  shown  by  figure 
2.5-   Figure  3.20  is  the  horizontal  cross  section  of  salinity 
anomaly,  (S7M-34.0%0  X  100,  at  a  depth  of  five  meters. 
Although  the  net  rate  of  evaporation  minus  precipitation 
curve  was  symmetric  about  the  equator,  the  surface  salinity 
pattern  was  asymmetric,  which  showed  the  strong  role  that 
advection  and  dynamics  played.   As  shown  by  figure  2.5>  two 
minima  of  evaporation  rate  minus  precipitation  rate  were 
located  at  about  6  N  and  6  S.   If  the  ocean  were  completely 
motionless,  two  minima  of  salinity  anomaly  would  tend  to 
result  at  these  locations,  as  the  net  rate  of  precipitation 
excess  would  tend  to  dilute  the  surface  water.   As  seen  from 
figure  3.20,  both  minima  were  shifted  poleward  by  the 
meridional  current  components  of  both  hemispheres.   Between 
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these  two  minima  a  relative  maximum  centered  at  approximately 
1  S  formed  due  to  the  relative  maximum  in  the  prescribed 
downward  flux  of  salt  and  a  contribution  due  to  vertical 
advection  at  the  eastern  boundary.   The  salinity  minimum  of 
the  Northern  Hemisphere  was  more  intense  than  its  southern 
hemispheric  counterpart,  a  feature  that  must  be  related  to 
the  asymmetric  surface  stress  X*,  (see  figure  2.3)  as  all 
other  boundary  conditions  were  symmetric  about  the  equator. 
Figure  3.21,  the  equatorial  E-W  salinity  anomaly  cross 
section,  further  showed  the  effect  of  salinity  advection. 
Since  the  prescribed  surface  flux  of  salt  in  the  equatorial 
regions  is  upward,  salinity  must  generally  increase  with 
depth.   Thus,  the  top  300  meters  are  generally  less  saline 
than  the  deeper  layers.   Within  this  upper  region  of  lower 
salinity  water  there  is  a  shallow  region  between  50  m  and 
100  m  of  even  less  saline  water  due  to  the  horizontal  conver- 
gence of  the  meridional  flux  of  fresh  water  from  north  and 
south  of  the  equator  (see  figure  3.23).   The  existence  of  a 
core  of  minimum  salinity  above  the  axis  of  the  countercurrent 
does  not  agree  with  the  results  of  Krivelevich  (1970)  in  his 
M0  layer  model.   He  found  a  core  of  minimum  salinity  in  the 
near  surface  region  and  reported  that  this  agreed  with 
observational  data. 

F.   DENSITY  STRUCTURE 

Figure  3-22  is  an  equatorial  E-W  cross  section  of  density 
anomaly  (£-  1.0)  K   1000.   Its  similarity  to  figure  3-18, 
the  corresponding  cross  section  for  temperature,  and  its 
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dissimilarity  to  figure  3.21,  especially  in  the  upper  100 
meters,  suggested  that  temperature  primarily  determined  the 
density  structure.   This  same  conclusion  was  reached  by 
Krivelevich  (1970)  who  used  a  nonlinear  equation  of  state. 
Because  of  the  strong  dependence  of  density  upon  temperature 
and  since  the  ocean  was  heated  from  above  which  tended  to 
insure  vertical  stability,  it  was  concluded  that  the  convec- 
tive  adjustment  process  played  only  a  minor  role  in  deter- 
mining the  density  structure.   Convective  adjustment  was 
thought  to  be  important  only  along  the  boundaries  where 
significant  vertical  velocities  existed. 

G.   MASS  TRANSPORT 

As  mentioned  previously,  if  the  currents  were  nearly 
geostrophic  in  the  vicinity  of  the  equator,  then  the  mass 
transport  across  the  equator  would  be  zero  at  each  level. 
Figure  3.23,  the  transport  streamfunction,  showed  however, 
there  was  transport  from  the  Northern  Hemisphere  to  the 
Southern  Hemisphere  at  the  surface  with  an  opposite  trans- 
port below.   The  dominant  flow  was  a  two  cell  circulation 
in  which  surface  water  was  transported  poleward  in  each 
hemisphere  where  it  struck  the  boundaries,  was  forced  to 
descend,  and  returned  as  deep  water  to  the  equator. 
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IV.   CONCLUSIONS 

The  model  is  capable  of  simulating  many  features  of  the 
observed  equatorial  ocean  flow  such  as  vertical  upwelling  at 
the  equator,  Ekman  divergence  and  surface  mass  transport 
towards  each  pole,  and  the  existence  of  a  subsurface  counter- 
current  with  a  realistic  meridional  profile.   The  failure  of 
the  model  to  accurately  resolve  the  countercurrent  profile 
below  100  m  was  probably  caused  by  using  a  vertical  eddy 
diffusion  coefficient  for  momentum,  k\*»    that  was  independent 
of  vertical  stability,  vertical  current  shear,  and  depth. 
A  better  modeling  of  vertical  mixing  by  using  a  variable  ^m 
will  probably  give  improved  results. 

The  salinity  pattern  was  shown  to  be  strongly  influenced 
by  the  model  currents  which  were  primarily  wind  driven, 
while  the  salinity  distribution  itself  appeared  to  have 
little,  if  any,  effect  on  the  general  current  circulations. 

The  failure  of  the  model  to  reach  temperature  or  salinity 
steady  state  for  the  bottom  layers  over  30  years  of  time 
integration  was  not  considered  serious  for  an  ocean  fore- 
casting model,  where  most  of  the  variation  in  parameters 
tends  to  occur  in  the  near  surface  layers.   However,  the 
slow  thermal  response  time  probably  helped  to  prevent  the 
formation  of  the  deep  isothermal  layer  and  some  means  of 
accurately  specifying  the  initial  temperature  profile  is 
clearly  needed. 
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